{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from osgeo import gdal, ogr, osr\n",
    "from pprint import pprint\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def raster2array(raster_fn):\n",
    "    raster = gdal.Open(raster_fn)\n",
    "    band = raster.GetRasterBand(1)\n",
    "    no_data = band.GetNoDataValue()\n",
    "    array = band.ReadAsArray()\n",
    "\n",
    "    return array, no_data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def point_index(array, no_data):\n",
    "    coordinates = []\n",
    "    for item in array[array != no_data]:\n",
    "        index = np.where(array == item)\n",
    "        coordinates.append(list(zip(index[0], index[1])))\n",
    "\n",
    "    return coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pixel2geo(raster_fn, point):\n",
    "    raster = gdal.Open(raster_fn)\n",
    "    band = raster.GetRasterBand(1)\n",
    "\n",
    "    geotransform = raster.GetGeoTransform()\n",
    "    originX = geotransform[0]\n",
    "    originY = geotransform[3]\n",
    "    pixelWidth = geotransform[1]\n",
    "    pixelHeight = geotransform[5]\n",
    "\n",
    "    x_offset, y_offset = point\n",
    "\n",
    "    x = x_offset * pixelWidth + originX\n",
    "    y = y_offset * pixelHeight + originY\n",
    "    return (x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_distance(point_a, point_b):\n",
    "    a_x, a_y = point_a\n",
    "    b_x, b_y = point_b\n",
    "\n",
    "    return math.sqrt(math.pow(b_x - a_x, 2) + math.pow(b_y - a_y, 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_minial_distance(target_point, coordinates):\n",
    "    p_x, p_y = target_point\n",
    "    \n",
    "    minimal = None\n",
    "    for point in coordinates:\n",
    "        distance = compute_distance(target_point, point)\n",
    "\n",
    "        if minimal is None:\n",
    "            minimal = distance\n",
    "        elif distance < minimal:\n",
    "            minimal = distance\n",
    "\n",
    "    return minimal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "def main(raster_fn):\n",
    "    array, no_data = raster2array('source_PointToRaster.tif')\n",
    "    index = point_index(array, no_data)\n",
    "\n",
    "    distance_array = np.zeros_like(array)\n",
    "\n",
    "    coordinates = []\n",
    "    for items in index:\n",
    "        for item in items:\n",
    "            point = pixel2geo(raster_fn, item)\n",
    "            coordinates.append(point)\n",
    "            \n",
    "            pprint(item)\n",
    "\n",
    "    for x, y in np.ndindex(array.shape):\n",
    "        geo_point = pixel2geo(raster_fn, (x, y))\n",
    "        distance_array[x, y] = find_minial_distance(geo_point, coordinates)\n",
    "    \n",
    "    return distance_array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0, 276)\n",
      "(35, 39)\n",
      "(78, 228)\n",
      "(82, 0)\n",
      "(114, 278)\n",
      "(133, 114)\n",
      "(206, 215)\n",
      "(207, 24)\n",
      "(282, 253)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[4265, 4204, 4145, ...,    0,   81,  162],\n",
       "       [4211, 4150, 4089, ...,   81,  115,  181],\n",
       "       [4158, 4096, 4035, ...,  162,  181,  230],\n",
       "       ...,\n",
       "       [6254, 6229, 6205, ..., 1879, 1960, 2041],\n",
       "       [6331, 6307, 6283, ..., 1873, 1955, 2036],\n",
       "       [6409, 6384, 6361, ..., 1871, 1953, 2034]], dtype=uint32)"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "main('source_PointToRaster.tif')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(283, 279)"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x, _ = raster2array('source_PointToRaster.tif')\n",
    "x.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "def offset2coordinate(geotransform, offset):\n",
    "    originX = geotransform[0]\n",
    "    originY = geotransform[3]\n",
    "    pixelWidth = geotransform[1]\n",
    "    pixelHeight = geotransform[5]\n",
    "\n",
    "    x_offset, y_offset = offset\n",
    "\n",
    "    x = x_offset * pixelWidth + originX\n",
    "    y = y_offset * pixelHeight + originY\n",
    "    return (x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "def array2raster_shp(shape_fn, geotransform, array):\n",
    "    shapefile = ogr.Open(shape_fn)\n",
    "    layer = shapefile.GetLayer()\n",
    "\n",
    "    rows, cols = array.shape\n",
    "\n",
    "    driver = gdal.GetDriverByName('GTiff')\n",
    "    out_raster = driver.Create('OUT.tif', cols, rows, 1, gdal.GDT_Float32)\n",
    "    out_raster.SetGeoTransform(geotransform)\n",
    "    out_band = out_raster.GetRasterBand(1)\n",
    "    out_band.WriteArray(array)\n",
    "\n",
    "    srs = layer.GetSpatialRef()\n",
    "    out_raster.SetProjection(srs.ExportToWkt())\n",
    "    out_band.FlushCache()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "maxdist = 1000\n",
    "distunits = 81.39"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "'304, 308'\n"
     ]
    }
   ],
   "source": [
    "data_source = ogr.Open('sourcen.shp')\n",
    "layer = data_source.GetLayer()\n",
    "\n",
    "point_x = []\n",
    "point_y = []\n",
    "points = []\n",
    "for feature in layer:\n",
    "    geom = feature.GetGeometryRef()\n",
    "    x, y, _ = geom.GetPoint()\n",
    "\n",
    "    point_x.append(x)\n",
    "    point_y.append(y)\n",
    "\n",
    "    points.append((x, y))\n",
    "\n",
    "x_min = min(point_x) - maxdist\n",
    "x_max = max(point_x) + maxdist\n",
    "y_min = min(point_y) - maxdist\n",
    "y_max = max(point_y) + maxdist\n",
    "\n",
    "rows = int((x_max - x_min) / distunits) + 1\n",
    "cols = int((y_max - y_min) / distunits) + 1\n",
    "\n",
    "pprint(f\"{rows}, {cols}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(304, 308)"
      ]
     },
     "execution_count": 58,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "distance_array = np.zeros((rows, cols))\n",
    "distance_array.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "array([[5686.64297919, 5631.27646667, 5576.54817783, ..., 7592.50748568,\n",
      "        7667.39741371, 7742.41854614],\n",
      "       [5627.55402255, 5571.60035789, 5516.28004035, ..., 7560.84889083,\n",
      "        7636.04931425, 7711.37542683],\n",
      "       [5569.02763908, 5512.47994247, 5456.56027656, ..., 7529.9369615 ,\n",
      "        7605.44303928, 7681.06930496],\n",
      "       ...,\n",
      "       [1390.63970587, 1333.31283788, 1278.59928044, ..., 2964.13028812,\n",
      "        2987.9179157 , 3013.7165164 ],\n",
      "       [1448.38167981, 1393.43238893, 1341.17367213, ..., 3042.35672746,\n",
      "        3065.53740202, 3090.68816503],\n",
      "       [1508.30661786, 1455.62123699, 1405.67641436, ..., 3120.74497733,\n",
      "        3143.34762945, 3167.88069375]])\n"
     ]
    }
   ],
   "source": [
    "# target_ds = gdal.GetDriverByName('MEM').Create('', rows, cols, 1, gdal.GDT_Float32)\n",
    "# target_ds.SetGeoTransform((x_min, distunits, 0, y_max, 0, -distunits))\n",
    "    \n",
    "for x, y in np.ndindex(distance_array.shape):\n",
    "    geo_point = offset2coordinate((x_min, distunits, 0, y_max, 0, -distunits), (x, y))\n",
    "    distance_array[x, y] = find_minial_distance(geo_point, points)\n",
    "\n",
    "for x, y in points:\n",
    "    x_offset = int((x - x_min) / distunits)\n",
    "    y_offset = int((y - y_max) / -distunits)\n",
    "    \n",
    "    distance_array[x_offset, y_offset] = 0\n",
    "    \n",
    "pprint(distance_array)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "array2raster_shp('sourcen.shp', pytho, distance_array.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
